Regression of Clinical Outcomes to Sigs

library(Seurat)
library(DAAG)
library(tidyverse)
library(relaimpo)
library(bootstrap)

Read in tumor object

subset tumor seurat obeject to TN only

tn_samples <- filter(tiss_subset_tumor2@meta.data, sample_name == "LT_S34" | sample_name == "LT_S43" | sample_name == "LT_S45" | sample_name == "LT_S49" | sample_name == "LT_S52" | sample_name == "LT_S51" | sample_name == "LT_S56" | sample_name == "LT_S67" | sample_name == "LT_S69" | sample_name == "LT_S74" | sample_name == "LT_S75")
tn_seurat <- SubsetData(tiss_subset_tumor2, cells.use = tn_samples$cell_id)
rownames(tn_seurat@meta.data) <- tn_seurat@meta.data$cell_id

Investigate each Signature found from grouped analysis: 1. Alveolar Sig 2. Kynurenine Sig 3. Plasminogen Sig 4. Serpine1 5. Gap Junction Sig

1. Alveolar Sig

DOR_Alveolar <- as.data.frame(FetchData(object = tn_seurat, vars.all = c("SFTPC", "SFTPB", "SFTPD", "PGC", "CLDN18", "AQP4", "SCGB3A1", "ABCA3", "GATA6", "NKX2-1", "SFTA3", "IGFBP2", "HOPX", "NAPSA", "FOXA2", "AGER", "LAMP1")))
DOR_Alveolar$cell_id <- rownames(DOR_Alveolar)
DOR_Alveolar <- merge(tn_seurat@meta.data, DOR_Alveolar, by = "cell_id")
rownames(DOR_Alveolar) <- DOR_Alveolar$cell_id

2. Kynurenine Sig

DOR_Kynurenine <- as.data.frame(FetchData(object = tn_seurat, vars.all = c('IDO1', 'KYNU', 'QPRT')))
DOR_Kynurenine$cell_id <- rownames(DOR_Kynurenine)
DOR_Kynurenine <- merge(tn_seurat@meta.data, DOR_Kynurenine, by = "cell_id")
rownames(DOR_Kynurenine) <- DOR_Kynurenine$cell_id

3. Plasminogen Sig

DOR_Plasminogen <- as.data.frame(FetchData(object = tn_seurat, vars.all = c('ANXA2', 'PLAT', 'PLAU', 'PLAUR')))
DOR_Plasminogen$cell_id <- rownames(DOR_Plasminogen)
DOR_Plasminogen <- merge(tn_seurat@meta.data, DOR_Plasminogen, by = "cell_id")
rownames(DOR_Plasminogen) <- DOR_Plasminogen$cell_id

4. Serpine1

DOR_SERPINE1 <- as.data.frame(FetchData(object = tn_seurat, vars.all = c('SERPINE1')))
DOR_SERPINE1$cell_id <- rownames(DOR_SERPINE1)
DOR_SERPINE1 <- merge(tn_seurat@meta.data, DOR_SERPINE1, by = "cell_id")
rownames(DOR_SERPINE1) <- DOR_SERPINE1$cell_id

5. Gap Junction Sig

DOR_GapJunction <- as.data.frame(FetchData(object = tn_seurat, vars.all = c('GJB3', 'GJB2', 'GJB4','GJB5')))
DOR_GapJunction$cell_id <- rownames(DOR_GapJunction)
DOR_GapJunction <- merge(tn_seurat@meta.data, DOR_GapJunction, by = "cell_id")
rownames(DOR_GapJunction) <- DOR_GapJunction$cell_id

fit 1 = Alveolar Sig

summary(fit1) # show results

Call:
lm(formula = dor ~ SFTPC + SFTPB + SFTPD + PGC + CLDN18 + AQP4 + 
    SCGB3A1 + ABCA3 + GATA6 + `NKX2-1` + SFTA3 + IGFBP2 + HOPX + 
    NAPSA + FOXA2 + AGER + LAMP1, data = DOR_Alveolar)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.31790 -0.03572  0.00897  0.04640  0.19908 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.6616254  0.0056816 116.451  < 2e-16 ***
SFTPC        0.0116742  0.0063559   1.837 0.066646 .  
SFTPB        0.0022941  0.0025315   0.906 0.365098    
SFTPD       -0.0429224  0.0107056  -4.009  6.7e-05 ***
PGC         -0.0082526  0.0135026  -0.611 0.541261    
CLDN18       0.0007957  0.0196329   0.041 0.967680    
AQP4        -0.0040422  0.0107509  -0.376 0.707031    
SCGB3A1     -0.0011107  0.0060576  -0.183 0.854568    
ABCA3       -0.0034290  0.0096255  -0.356 0.721763    
GATA6       -0.0226240  0.0228514  -0.990 0.322471    
`NKX2-1`    -0.0451323  0.0044903 -10.051  < 2e-16 ***
SFTA3        0.0054035  0.0059639   0.906 0.365205    
IGFBP2       0.0330716  0.0016918  19.548  < 2e-16 ***
HOPX        -0.0335134  0.0036519  -9.177  < 2e-16 ***
NAPSA       -0.0302077  0.0033578  -8.996  < 2e-16 ***
FOXA2       -0.0273386  0.0044887  -6.091  1.8e-09 ***
AGER        -0.0907966  0.0243912  -3.723 0.000212 ***
LAMP1        0.0030322  0.0095368   0.318 0.750610    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.071 on 745 degrees of freedom
Multiple R-squared:  0.7277,    Adjusted R-squared:  0.7214 
F-statistic: 117.1 on 17 and 745 DF,  p-value: < 2.2e-16

fit2 = Kynurenine Sig

fit2 <- lm(dor ~ IDO1 + KYNU + QPRT, data=DOR_Kynurenine)
summary(fit2) # show results

Call:
lm(formula = dor ~ IDO1 + KYNU + QPRT, data = DOR_Kynurenine)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.35619 -0.09619 -0.09157  0.15381  0.25432 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.656187   0.005075 129.289  < 2e-16 ***
IDO1        -0.025909   0.013803  -1.877  0.06089 .  
KYNU         0.026314   0.009387   2.803  0.00519 ** 
QPRT        -0.103465   0.015602  -6.631 6.31e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1303 on 759 degrees of freedom
Multiple R-squared:  0.0659,    Adjusted R-squared:  0.0622 
F-statistic: 17.85 on 3 and 759 DF,  p-value: 3.348e-11

# diagnostic plots 
plot(fit2)

ggplot(DOR_Kynurenine, aes(x = IDO1, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_Kynurenine, aes(x = KYNU, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_Kynurenine, aes(x = QPRT, y = dor, color = sample_name)) + geom_point()

fit3 = Plasminogen Sig

fit3 <- lm(dor ~ PLAU + PLAUR + PLAT + ANXA2, data=DOR_Plasminogen)
summary(fit3) # show results

Call:
lm(formula = dor ~ PLAU + PLAUR + PLAT + ANXA2, data = DOR_Plasminogen)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.39284 -0.07108  0.02876  0.07998  0.24803 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.515175   0.013690  37.631  < 2e-16 ***
PLAU        -0.026385   0.004716  -5.595 3.08e-08 ***
PLAUR       -0.031850   0.006093  -5.227 2.22e-07 ***
PLAT        -0.043839   0.003904 -11.229  < 2e-16 ***
ANXA2        0.056471   0.003923  14.393  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1065 on 758 degrees of freedom
Multiple R-squared:  0.3766,    Adjusted R-squared:  0.3733 
F-statistic: 114.5 on 4 and 758 DF,  p-value: < 2.2e-16

# diagnostic plots 
plot(fit3)

ggplot(DOR_Plasminogen, aes(x = PLAU, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_Plasminogen, aes(x = PLAUR, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_Plasminogen, aes(x = PLAT, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_Plasminogen, aes(x = ANXA2, y = dor, color = sample_name)) + geom_point()

fit4 = SERPINE1

fit4 <- lm(dor ~ SERPINE1, data=DOR_SERPINE1)
summary(fit4) # show results

Call:
lm(formula = dor ~ SERPINE1, data = DOR_SERPINE1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.36405 -0.10405 -0.04199  0.14595  0.25223 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.664046   0.005113 129.880  < 2e-16 ***
SERPINE1    -0.047593   0.007384  -6.446 2.04e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1311 on 761 degrees of freedom
Multiple R-squared:  0.05177,   Adjusted R-squared:  0.05052 
F-statistic: 41.55 on 1 and 761 DF,  p-value: 2.042e-10

# diagnostic plots 
plot(fit4)

ggplot(DOR_SERPINE1, aes(x = SERPINE1, y = dor, color = sample_name)) + geom_point()

fit5 = Gap Junction Sig

fit5 <- lm(dor ~ GJB3 + GJB2 + GJB4 + GJB5, data=DOR_GapJunction)
summary(fit5) # show results

Call:
lm(formula = dor ~ GJB3 + GJB2 + GJB4 + GJB5, data = DOR_GapJunction)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.33107 -0.07003 -0.07003  0.14248  0.17997 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.630034   0.004815 130.836  < 2e-16 ***
GJB3        0.029685   0.018856   1.574   0.1158    
GJB2        0.042758   0.023258   1.838   0.0664 .  
GJB4        0.102237   0.016175   6.321 4.44e-10 ***
GJB5        0.124777   0.014473   8.622  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1227 on 758 degrees of freedom
Multiple R-squared:  0.1722,    Adjusted R-squared:  0.1678 
F-statistic: 39.42 on 4 and 758 DF,  p-value: < 2.2e-16

# diagnostic plots 
plot(fit5)

ggplot(DOR_GapJunction, aes(x = GJB2, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_GapJunction, aes(x = GJB3, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_GapJunction, aes(x = GJB4, y = dor, color = sample_name)) + geom_point()

ggplot(DOR_GapJunction, aes(x = GJB5, y = dor, color = sample_name)) + geom_point()

# # K-fold cross-validation
# cv.lm(data = DOR_GapJunction, form.lm = fit5, m = 10, plotit = FALSE)
# # Assessing R2 shrinkage using 10-Fold Cross-Validation 
# # define functions 
# theta.fit <- function(x,y){lsfit(x,y)}
# theta.predict <- function(fit5,x){cbind(1,x)%*%fit5$coef} 
# 
# # matrix of predictors
# X <- as.matrix(DOR_GapJunction[c("GJB3","GJB2","GJB4","GJB5")])
# # vector of predicted values
# y <- as.matrix(DOR_GapJunction[c("dor")]) 
# 
# results <- crossval(X,y,theta.fit,theta.predict,ngroup=10)
# cor(y, fit5$fitted.values)**2 # raw R2 
# cor(y,results$cv.fit5)**2 # cross-validated R2
# 
# # Calculate Relative Importance for Each Predictor
# calc.relimp(fit5,type = c("lmg","last","first","pratt"), rela=TRUE)
# # Bootstrap Measures of Relative Importance (1000 samples) 
# boot <- boot.relimp(fit5, b = 1000, type = c("lmg", "last", "first", "pratt"), rank = TRUE, diff = TRUE, rela = TRUE)
# booteval.relimp(boot) # print result
# plot(booteval.relimp(boot,sort=TRUE)) # plot result

table(tn_seurat@meta.data$biopsy_site, tn_seurat@meta.data$dor)

         
          0.3 0.31 0.43 0.46 0.5 0.56 0.57 0.7 0.81
  Adrenal   0    1    0    0   0    0    0   0    0
  Brain     0    0    0    0   0    0    0   0    0
  Liver     0    0    0    0   0    0    0  28    0
  LN        0    0    5    0   6  305   16   0    0
  Lung     14    0    0    0   0   71    0   0  293
  Pleura    0    0    0   24   0    0    0   0    0

table(tn_seurat@meta.data$sample_name, tn_seurat@meta.data$dor)

        
         0.3 0.31 0.43 0.46 0.5 0.56 0.57 0.7 0.81
  LT_S34   0    0    0   24   0    0    0   0    0
  LT_S43   0    0    5    0   0    0    0   0    0
  LT_S45   0    1    0    0   0    0    0   0    0
  LT_S49   0    0    0    0   6    0    0   0    0
  LT_S51   0    0    0    0   0    0   16   0    0
  LT_S52  14    0    0    0   0    0    0   0    0
  LT_S56   0    0    0    0   0    0    0   0  291
  LT_S67   0    0    0    0   0    0    0   0    2
  LT_S69   0    0    0    0   0  305    0   0    0
  LT_S74   0    0    0    0   0   71    0   0    0
  LT_S75   0    0    0    0   0    0    0  28    0

table(tn_seurat@meta.data$sample_name)

LT_S34 LT_S43 LT_S45 LT_S49 LT_S51 LT_S52 LT_S56 LT_S67 LT_S69 LT_S74 LT_S75 
    24      5      1      6     16     14    291      2    305     71     28 

Bulkize the samples

tn_seurat <- SetIdent(tn_seurat, ident.use = tn_seurat@meta.data$sample_name)
table(tn_seurat@ident)

LT_S34 LT_S43 LT_S45 LT_S49 LT_S51 LT_S52 LT_S56 LT_S67 LT_S69 LT_S74 LT_S75 
    24      5      1      6     16     14    291      2    305     71     28 

sample.averages <- AverageExpression(object = tn_seurat)

Finished averaging RNA for cluster LT_S34
Finished averaging RNA for cluster LT_S43
Finished averaging RNA for cluster LT_S45
Finished averaging RNA for cluster LT_S49
Finished averaging RNA for cluster LT_S51
Finished averaging RNA for cluster LT_S52
Finished averaging RNA for cluster LT_S56
Finished averaging RNA for cluster LT_S67
Finished averaging RNA for cluster LT_S69
Finished averaging RNA for cluster LT_S74
Finished averaging RNA for cluster LT_S75

To find DE genes between bulkized TN samples with low and high DOR, export table with groups

# set up table 
sample.averages.t <- as.data.frame(t(sample.averages))
head(sample.averages.t)

sample.averages.t$sample_name <- rownames(sample.averages.t)
sample.averages.t <- left_join(sample.averages.t, dor_meta, by = "sample_name")
rownames(sample.averages.t) <- sample.averages.t$sample_name
length(colnames(sample.averages.t))

[1] 26489

DE_avg <- pairwise.wilcox.test(x = sample.averages.t$EGFR, g = sample.averages.t$dor_class)
write.csv(sample.averages.t, file = "/myVolume/TN_bulkized_data.csv")
TN.sample.averages <- sample.averages
head(TN.sample.averages)

Bulkize fit analysis Alveolar

Bulkize fit analysis Kynurenine

Bulkize fit analysis Plasminogen

Bulkize fit analysis of SERPINE1

Bulkize fit analysis GapJunction

bulkized_TN_markers <- read.csv(file = paste(dir, "Data_input/mwu_luad.csv", sep = ""))
bulkized_TN_markers.f <- filter(bulkized_TN_markers, pval_1 <= 0.05)
hist(bulkized_TN_markers.f$stat_1)

length(bulkized_TN_markers.f$pval_1)

[1] 4115

bulkized_TN_markers.f <- bulkized_TN_markers.f[order(bulkized_TN_markers.f$stat_1, decreasing = FALSE), ] 

Most compelling high expression corr to low dor

ggplot(sample.averages.t, aes(x = ADAR, y = dor)) + geom_point(aes(color = patient_id))

ggplot(sample.averages.t, aes(x = CFL1, y = dor)) + geom_point(aes(color = patient_id))

Most compelling high expression corr to high dor

ggplot(sample.averages.t, aes(x = TTLL13P, y = dor)) + geom_point(aes(color = dor_class))

ggplot(sample.averages.t, aes(x = ALS2, y = dor)) + geom_point(aes(color = dor_class))

ggplot(sample.averages.t, aes(x = RLN1, y = dor)) + geom_point(aes(color = dor_class))

ggplot(sample.averages.t, aes(x = USP45, y = dor)) + geom_point(aes(color = dor_class))

ggplot(sample.averages.t, aes(x = BDKRB1, y = dor)) + geom_point(aes(color = dor_class))

ggplot(sample.averages.t, aes(x = LINC01061, y = dor)) + geom_point(aes(color = dor_class))

ggplot(sample.averages.t, aes(x = ZNF563, y = dor)) + geom_point(aes(color = dor_class))

ggplot(sample.averages.t, aes(x = WDR19, y = dor)) + geom_point(aes(color = dor_class))